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Abstract. We investigate the dynamics and evolution of coa- 
scing neutron stars. The three-dimensional Newtonian equa- 
ins of hydrodynamics are integrated by the "Piecewise 
pbohc Method" on an equidistant Cartesian grid with a res- 
ation of 64-' or 128^. Although the code is purely Newtonian, 
do include the emission of gravitational waves and their 
^j^kreaction on the hydrodynamic flow. The properties of neu- 
trc^n star matter are described by the physical equation of state 
of Lattimer & Swesty(1991).In addition to the fundamental hy- 
tpdynamic quantities, density, momentum, and energy, we fol- 
' the time evolution of the electron density in the stellar gas. 
jergy loss by all types of neutrinos and changes of the electron 
^J^ction due to the emission of electron neutrinos and antineu- 
^^os are taken into account by an elaborate "neutrino leakage 
Jscneme". We simulate the coalescence of two identical, cool 
^eutron stars with a baryonic mass of « 1 .6 Mq and a radius of 
rSis km and with an initial center-to-center distance of 42 km. 
^fe initial distributions of density and electron concentration 
^ given from a model of a cold neutron star in hydrostatic 
^^i^aiUbrium, the temperature in our initial models is increased 

th that the thermal energy is about 3% of the degeneracy 
rgy for given density and electron fraction (central temper- 
ature about 8 MeV). We investigate three cases which differ by 
the initial velocity distribution in the neutron stars, representing 
different cases of the neutron star spins relative to the direction 
of the orbital angular momentum vector. The orbit decays due to 
gravitational-wave emission, and after half a revolution the stars 
are so close that dynamical instability sets in. Within about 1 ms 
they merge into a rapidly spinning (Pspin « 1 ms), high-density 
body (p « 10'^ g/cnr') with a surrounding thick disk of mate- 
rial with densities p « 10^*^ — lO'^ g/cm^ and orbital velocities 
of 0.3-0.5 c. In this work we evaluate the models in detail with 
respect to the gravitational wave enussion using the quadrupole 



approximation. In a forthcoming paper we wiU concentrate on 
the neutrino emission and implications for gannma-ray bursters. 
The peak emission of gravitational waves is short but powerful. 
A maximum luminosity in excess of 10^^ erg/s is reached for 
about 1 ms. The amphtudes of the gravitational waves are close 
to 3 • 10~^^ at a distance of 1 Gpc, and the typical frequencies 
are between 1 KHz and 2 KHz, near the dynamical frequency of 
the orbital motion of the merging and coalescing neutron stars. 
In contrast to the diverging gravitational wave ampUtude of two 
coalescing point-masses, our models show decreasing amph- 
tudes of the waves because of the finite extension of the neutron 
stars and the nearly spherical shape of the merged object toward 
the end of the simulations. The structure and temporal devel- 
opment of the gravitational wave signal and energy spectrum 
show systematic trends with the amount of angular momentum 
in the system and depend on the details of the hydrodynamic 
mass motions. 
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1. Introduction 

Gravitational wave astronomy may soon become an obser- 
vational science, since three gravitational wave experiments, 
LIGO (Abramovici et al, 1992), VIRGO (Bradaschia et al, 
1991), and GEO600 (Danzmann et al, 1995a) with four in- 
terferometric detectors are under construction. Moreover, the 
European Space Agency (ESA) has chosen the space-based 
LISA laser interferometer (Danzmann et al, 1995b) to be a cor- 
nerstone project. All sources are expected to be astrophysical, 
with merging and coalescence of neutron stars and black holes 
as well as rapidly rotating neutron stars and supemovae being 
the most promising candidates (Thorne, 1992). Whether these 
events and objects will be detectable or not does not only depend 
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on their occurrence rates and the strength of the emitted waves, 
but it also depends on the frequency of the gravitational waves 
relative to the resonance frequencies of the detectors (e.g. / « 
10-1000 Hz for terrestrial detectors) (Finn & Chernoff, 1993; 
Finn, 1994). Only an a priori knowledge of the structure of 
typical gravitational wave signals will allow one to extract the 
faint signals from a very noisy background. The event rates 
are highly uncertain and depend on the Hubble constant, pulsar 
and supernova rates, and orbital parameters of binary systems. 
Very optimistic estimates yield several neutron star and black 
hole mergers per year within about 25 Mpc, while pessimistic 
numbers quote a few per year within 1 Gpc (e.g. Clark et al, 
1979; Phinney, 1991). 

Two directions of theoretical investigations have been fol- 
lowed up. (a) Calculations of the orbital decay of binaries due to 
gravitational radiation, treating the binary components as point- 
masses and using progressively more elaborate post-Newtonian 
approximations to the equations of general relativity. These 
computations can integrate the orbit over lots of revolutions. Ex- 
amples of such investigations are, e.g., Lincoln & Will (1990), 
Iyer & Will, (1993), Cutler & Hanagan, (1994), Imshennik & 
Popov (1994). (b) Hydrodynamic simulations that concentrate 
on the merging phase of the binary and the last few orbital 
revolutions before the merging. These simulations take into ac- 
count the finite extension of the binary components and allow 
one to study the fluid dynamics effects and to include detailed 
microphysics. In particular, computations with different numer- 
ical methods, with different handling of the gravitational wave 
emission and of the backreaction on the hydrodynamic flow, 
and with different parameters in a polytropic equation of state 
have been performed so far (e.g., Shibata, Nakamura & Oohara, 
1993; Zhuge, CentreUa & McMiUan 1994; Rasio & Shapiro 
1994; Davies et al, 1994; also references cited in these papers). 
Our work intends to go a step further in the numerical modeling 
of the hydrodynamic phase of the inspiral by using a complex, 
"physical" equation of state (Lattimer & Swesty, 1991) and by 
including the effects of local composition changes and cooling 
due to neutrino emission. 

Gravitational waveforms and luminosities were calculated 
in the post-Newtonian expansion for the small parameter 
e = (v Icf to terms of order 0{e^l^) by Lincoln & WiU (1990). 
Iyer & Will (1993) derived the post-Newtonian radiation reac- 
tion terms at Oie' 1^). Imshennik & Popov (1994) showed that 
independent of the initial parameters of the binary orbit the final 
orbital eccentricity gets very small. Cutler & Flanagan (1994) 
investigated which, and how accurately, binary parameters can 
be deduced from the waveforms. They found, e.g., that although 
to lowest post-Newtonian order the phase depends only on the 
"chirp mass' ' M = {M\M2f/^{Mi + Mz)-^/^ which IS mea- 
surable to a relative accuracy better than 1%, higher terms also 
involve the reduced mass. 

In a sequence of papers a Japanese group (Shibata, Naka- 
mura & Oohara, e.g., 1993) reported about hydrodynamic sim- 
ulations with different initial conditions of the neutron stars' 
spins, masses, and tidal deformations and with different shapes 
of the orbit. The maximum amplitude of the gravitational radia- 



tion at 10 Mpc turned out to be 1-510^ and the total amount 
of energy emitted in gravitational waves was determined to be 
1-3% of the total rest mass. Especially the latter values are 
strongly dependent on the size and compactness of the neu- 
tron stars since the energy in gravitational waves increases like 
1 /r^ with the minimum separation of the neutron stars before 
they merge and the binary system loses its highly non-spherical 
geometry. 

Zhuge et al (1994) focussed on the energy spectrum that 
gives the energy emitted in gravitational waves at different fre- 
quencies. They discussed the information about the dynamics 
of the coalescence and about the merging objects that may be 
provided by the spectrum. Certain stages of the binary evolution 
are reflected in particular wave features. From their analysis it 
is evident that the equation of state and the neutron star masses 
and radii influence the frequencies and amplitudes of spectral 
structures. Therefore an interpretation of details and fine struc- 
ture in gravitational wave signals seems to require models of 
the merger scenario that include a physically more meaningful 
description of neutron star matter than by the usually employed 
polytropic equations of state. Microphysical processes at high 
densities, like changes of the /3-equilibrium, give rise to bulk 
viscosity and neutrino emission and must be suspected to cause 
dissipation and to smooth out or modify some of the fine struc- 
ture found in idealized models. The same effect may be associ- 
ated with the gravitational radiation (back)reaction which was 
also not taken into account by Zhuge et al (1994). 

Rasio & Shapiro (1994) found that triaxial configurations 
are formed by the merging of neutron stars with a sufficiently 
stiff equation of state. Non-axisymmetric structures occur for 
adiabatic exponents of 3 or higher (polytropic index 0.5 or be- 
low) in the polytropic equation of state. In this case the peak am- 
pUtude of the gravitational wave emission is substantially larger 
and the emission proceeds for a longer time after the coales- 
cence. Using hydrodynamic models that were computed with a 
polytropic equation of state Davies et al (1994) tried to estimate 
the temperatures in the central, merged object and in the sur- 
rounding thick disk and obtained values of around 10 MeV and 
3 MeV, respectively. Based on these estimates they attempted 
to discuss possible implications for neutrino emission, gamma- 
ray bursts, and thermonuclear reactions and energy generation 
in the disk. 

The hydrodynamic simulations presented in this work em- 
ploy the high-density equation of state of Lattimer & Swesty 
(1991) which describes the thermodynamics of the neutron star 
matter in dependence of density, temperature, and composition, 
the latter being expressed by the electron fraction Y^. Changes 
of Ye and cooUng of the matter due to the emission of neutrinos 
of all types are taken into account. The use of a microphysically 
meaningful equation of state and the consideration of neutrino 
source terms imply that dissipation effects due to the bulk vis- 
cosity of the medium are included in our models. Although 
our computations, like those of the hydrodynamic approaches 
summarized above, were performed with a basically Newto- 
nian code, we add the terms into the hydrodynamic equations 
that describe the effects of gravitational wave emission and the 
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corresponding backreaction on the hydrodynamic flow. This is 
done according to the formalism of Blanche! et al (1990) by 
following the special implementation chosen by Shibata et al 
(1992). Certainly, not taking into account the gravitational ra- 
diation reaction (Rasio & Shapiro, 1994) is a poor choice. Our 
treatment is also more accurate than the point-mass quadrupole 
approximation that was used by Davies et al (1994) and Zhuge 
et al (1994) for the phase of the coalescence when the neutron 
stars were still separated and switched off when the merging 
happened. 

Our simulations are performed with an explicit, Eulerian 
finite-difference scheme based on the PPM-method of Colella 
& Woodward (1984). This grid-based scheme is superior to 
a particle-based method hke SPH in the handling of shock 
waves, which is done very well by solving local Riemann 
problems. Compared with other Eulerian algorithms the PPM 
method has also a rather small numerical viscosity. This seems 
to be demanded since recent work (Kochanek, 1992; Lai, 1994; 
Reisenegger & Goldreich, 1994; see Sect. 3 for a discussion) 
suggests that the viscosity of neutron star matter is so small that 
the stars cannot develop synchronous rotation during inspiral. 

Our models yield detailed information about the neutrino 
emission (lepton number as well as energy) from the merging 
stars. This allows us to analyse the merger scenario not only as 
a source of gravitational radiation but also to perform a quanti- 
tative investigation of the possibihty to power gamma-ray burst 
events by the energy deposition due to neutrino-antineutrino an- 
nihilation in the surroundings of the merging binary. We shall 
report results of this part of our work in a subsequent paper, 
while we will concentrate on the gravitational wave aspect here. 
Section 2 provides a description of the numerical methods and 
of the treatment of the microphysics. In Sect. 3 we explain the 
set of computed models and give details about the initial condi- 
tions our simulations are started from. The results are presented 
in Sect. 4, followed by a discussion in Sect. 5. 



2. Computational procedure 

2.1. Hydrodynamics, self- gravity and gravitational-wave back- 
reaction 

The scenario of coalescing neutron stars is a three-dimensional 
problem with the orbital plane being a plane of symmetry. We 
perform the simulations on an equidistant Cartesian grid with 
the number of zones being 64 or 128 in both space dimensions 
of the orbital plane. Perpendicular to the orbital plane (i.e. in 
the z-direction) we use only 1/4 of this number of zones. One 
factor of two is saved due to the symmetry about the orbital 
plane: we only simulate the volume on one side of this plane. 
A second factor of two can be saved by only using a region 
vertical to the orbital plane which has a quarter of the size 
of the region modeled in the orbital plane. This is possible 
because test calculations showed that hardly any matter moves 
out to more than one neutron star radius away from the orbital 
plane. 



The neutron star matter is evolved hydrodynamicly via the 
Piecewise Parabolic Method (PPM) developed by Colella and 
Woodward (1984). Details of the actual implementation of the 
PPM code can be found in Fryxell et al (1989) and in Ruf- 
fert (1992). It is Eulerian and explicit in time and includes 
the non-constant-7 equation of state formaUsm as described by 
Colella & Glaz (1985). The code is basically Newtonian, but 
was extended to contain the terms necessary to describe grav- 
itational waves and their backreaction (Blanchet et al, 1990). 
We use a form similar to Shibata et al (1992) in notation and 
in use of the energy equation instead of the enthalpy equa- 
tion. However, we employ the total specific energy, i.e. the sum 
of specific internal energy and specific kinetic energy as fun- 
damental hydrodynamic variable. The hydrodynamic laws of 
conservation of mass, momentum and energy for non- viscous 
flow, in ADM coordinates (generaUzed isotropic coordinates), 
are, respectively. 



dp dpv^ 

-TT + -^-^ = , 
dt dxi 

dpw^ dipw^v^ + PS^^) _ dtp d(j) 

dt dx^ ^ dx'^ ^ dx'^ 

dpE d{pE + P)vi J d,l; 

-m^ dx^ --p^^a^ + ^ + ^B 



(1) 

(2) 
(3) 



with the energy source term due to gravitational waves given 
by 

The symbols have the usual meanings: p denotes mass density, 
t time, 11* and the components of the kinematic and dy- 
namic velocity, respectively, position vector components, P 
pressure, the Newtonian gravitational potential, cf) the backre- 
action potential due to gravitational waves, E is the total specific 
energy (sum of specific internal energy e and specific kinetic en- 
ergy jw^w^), G and c the gravitational constant and the speed 
of light, Djj the third time derivative of the quadrupole mo- 
ment, and 5e the energy loss due to neutrino emission. Details 
of S-E follow in Sect. 2.4 and Appendix B. The total emitted 
luminosity of gravitational waves C can be obtained either by 
summing up W over the whole emitting volume or via the 
classical quadrupole formula 



1 G 



(5) 



Notice that summing up W yields the gravitational wave lumi- 
nosity without averaging over time. When the orbit decays and 
therefore is not perfectly circular, the non-averaged luminosity 
is in general not identical with the value obtained by averaging 
over one orbital period. All secondary quantities are calculated 
from the primary quantities p, pw\ and pE (used in the code) 
by the following relations 



E = e + l-w^w'' 
2 



(6) 
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4G ... 



2 G 

= - ^ \R- Dij X- 



(7) 
(8) 



EOS, Ye=0.046 



Di 



STF 



I 



dV 



„ „ dvi dtp , 
ox^ ox^ 



dpv'' 
dx'' 



2pVi) - PXig — 



The notation STF means symmetric and frace free, 
111 



(9) 



(10) 



and the following three Poisson equations have to be solved: 



A4> = AttGp , 



Alp = -A-kG 



dx'' 



AR = 4tvG Dii x^ 



dx' 



(11) 



(12) 



(13) 



The Poisson equations in integral form are interpreted as con- 
volution and calculated by fast Fourier transform routines: 
non-periodic grid boundaries are enforced by zero-padding 
(e.g. Press et al, 1992; Eastwood & Brownrigg, 1979). This 
doubles the number of Fourier components compared to the 
original number of grid zones per dimension. The accelera- 
tions that follow from the potential are added as source terms 
in the PPM algorithm. All spatial derivatives necessary for the 
gravitational-wave terms are implemented as standard centered 
differences on the grid. 

We implement the hydrodynamic laws including gravita- 
tional backreaction (Eqs. 1 to 3) in the following way. First the 
PPM routines produce the kinematic velocities at the zone in- 
terfaces by solving Newtonian Riemann problems. Then these 
kinematic velocities are transformed to dynamic velocities 
(Eq. 7) and the relativistic momentum (pw*) and total specific 
energy E are advected. Finally the new relativistic quantities 
are transformed back to kinematic velocities to be used as start 
values for the next explicit PPM step. 

2.2. Equation of state 

In general, the state of matter in equilibrium is characterized 
by three independent quantities, e.g. baryon (rest mass) density 
p, electron fraction Yg, and internal energy density e = ep. All 
other thermodynamic quantities, e.g. pressure P, temperature 
T, etc., can be computed from these fundamental quantities. 
Since we include neutrino processes, the electron number den- 
sity He has to be advected, too: 



(14) 




20 30 
temperature T [MeV] 
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Fig. 1. An excerpt of the equation of state: the internal energy density 
as a function of temperature for three densities at one specific value of 



The electron number fraction Yg can be recovered from 
Ye = — u , 



(15) 



where u denotes the atomic mass unit. Information about the 
lepton source term Sl can be found in Sect. 2.4 and Ap- 
pendix B. 

We use the equation of state of Lattimer & Swesty (1991) 
in tabular form (incompressibility parameter of bulk nuclear 
matter K= 180 MeV). The tabulated quantities are the pres- 
sure P, the internal energy density e, the adiabatic index 
F = (dlnP/dlnp)\s, and the degeneracy parameters of pro- 
tons and of neutrons (without rest masses). The tables span the 
following ranges: (a) density: 9 ^ Igp [g/cm^] <> 15.5 in 130 
steps, (b) the temperature: -1.0 £ IgT [MeV] <, 1.7 in 108 
steps, and (c) the electron number fraction: 0.006 ^Yg^ 0.49 
in 25 steps. First a large table was generated with double as 
many entries in every dimension. The entries at small values 
of Ye were obtained from extrapolation, since the original Lat- 
timer & Swesty (1991) FORTRAN version is not appUcable 
for Ye ^ 0.03. Although Lattimer & Swesty (1991) modeled 
the equation of state to arbitrarily high densities, they did not 
consider possible additional physics such as pion or kaon con- 
densation, hyperons, and the quark-hadron phase transition. 
Therefore their equation of state provides a physical descrip- 
tion of nuclear matter only below about twice nuclear den- 
sity. Note, however, that only relatively small regions in our 
M « 1 .6 Mq neutron star models and in our simulated merg- 
ers actually reach densities above this value. If the FORTRAN 
version of the equation of state of Lattimer & Swesty (1991) 
did not yield converged values, the tabulated quantities were 
interpolated from neighbouring entries. Non-monotonic values 
of the energy dependence on temperature were smoothed by 
averaging over neighbouring entries in the table. Unambiguous 
inversion of the equation of state requires strictly monotonic 
behaviour Also isolated positive and negative spikes of the adi- 
abatic index F were reduced: negative ones by replacing the 
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entries with F = 0.01 and spikes over F = 4 were replaced 
by the average of six neighbouring entries in the table. Finally, 
the reduced table version (with the binning (a), (b), and (c) as 
specified above) was obtained by choosing new grid points in 
the Igp-lgT-Fe-space which are cubic-centered between the 
points of the large table and attributing to them the averages 
of the eight (=2-') closest neighbours of the large table. This 
reduces the amount of data and smooths the functions further. 

The basic quantities used by the PPM code are density 
p, electron concentration Y^, and energy E. Thus the tables 
have to be inverted to obtain the temperature. This is done 
by a bisection iteration. The interpolations in the tables are 
performed trilinearly in Ig p, Ig T and ■ A particular problem 
with the equation of state arises in the very degenerate, low- 
temperature regime, where the internal energy density e varies 
only weakly with temperature. Therefore small variations of e 
imply large shifts of T. Exemplary curves of e(T) are plotted 
in Fig. 1. 

The inversion of the table for typical values of density, en- 
ergy and Ye was checked by first computing the corresponding 
temperatures and then using the table a second time in reverse 
direction to regain the original quantities. The relative difference 
between these values and the original ones amounts to approxi- 
mately machine accuracy (10^^^) in all tested cases. Because of 
the numerical problems and the indispensable smoothing steps 
described above we consider the employed tabular resolution 
to be sufficient and to ensure the thermodynamical consistency 
of the equation of state to an acceptable degree. 

2.3. Gravitational wave emission 

Apart from the gravitational wave luminosity (Eq. 5), the am- 
plitudes and waveforms of the gravitational wave emission are 
of immediate interest, because they are the quantities that de- 
termine the measurements. The amplitudes of the two physical 
(i.e. coordinate-invariant or gauge-invariant) polarizations of 
gravitational waves that are emitted perpendicular to the orbital 
plane and are observed at a distance r are given by 



I, G \ (■■ ■■ \ 



and 



G2 



D 



xy 



(16) 



(17) 



In the direction along the rotation axis the gravitational-wave 
signal is largest. The amplitude h+ of the radiation in direction 
{6, if) then follows to (Rasio & Shapiro, 1994) 

h+{e,ip) = -h+iO,0) cos9 sinlip + hxiO,0) cos9 coslip , (18) 

where /i+(0,0) and /ix(0,0) are the amplitudes along the ro- 
tation axes, as given by Eqs. 16 and 17, respectively. hxiO, ip) 
cannot be simply written as function of h+{0, 0) and /ix (0, 0) 
only. 



We calculate the second time derivative of the quadrupole 
moment from 



Dij= STF 



(19) 



Neither in these equations, nor in those pertaining to the 
source terms necessary in the hydrodynamic code (Eqs. 4, 9, 
etc.), need time derivatives be calculated explicitly. Therefore, 
smoothing of gravitational wave quantities is unnecessary and 
all plots containing gravitational wave properties are done with- 
out any smoothing, except where noted. 

Analytic expressions for quasi-circular inspiral can be de- 
rived in the point-mass approximation (e.g. Misner et al, 1973). 
The separation a(t) of the two point masses, each with mass M, 
as a function of time t is 

/ / \ 1/4 5 4 

a{t) = ao ( 1 - - ) with to=—-^ . (20) 



6ARlc 



t() is a function of the separation ao at time i = 0. i?s = 2GM/c^ 
is the Schwarzschild-radius of mass M. The gravitational waves 
emitted then follow to (e.g. Lai et al, 1994) 



hli0 = O,ip,t): 



and 



hli0 = O,ip,t) = 



IK 

r ait) 



sin($(i) - (p) 



1 

raii) 



cos($(i) - (p) 



(21) 



(22) 



where the superscript p indicates the point-mass result. The 
angle $(i) as a function of time is given by 



*(i) = $0 



1 /a(i)\ 



5/2 



(23) 



where $0 = ^(t = to) is the final angle between the axis con- 
necting the two point masses and the direction ip = 0. The 
power radiated away in gravitational waves is (e.g. Shapiro & 
Teukolsky, 1983) 



5G 



(24) 



The gravitational wave luminosity of Eq. 24 represents a time- 
averaged value. Note that for circular orbits, the local power ex- 
pression Eq. 5 yields Eq. 24 even without averaging over one 
orbital period. Our results concerning the gravitational wave 
luminosity (Fig. 22) are somewhat coordinate dependent, espe- 
cially during the dynamic "plunge", although the total amount 
of gravitational wave energy emitted {£ in Table 1) as well as 
the energy spectrum (Fig. 27) are independent of the coordinate 
system. 

The energy emitted in gravitational waves per unit fre- 
quency interval d£/df is given (Zhuge et al, 1994; Thome, 
1987; Misner et al, 1973, Box 37.4.C) by 



- TT / 



d£ 

df c 



G 8 



5 15 

- C, 



(25) 



Czz I + \Cx 
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where C,j contains the Fourier transform of Dij. We numer- 
ically calculate the six power spectra via the discrete fast 
Fourier transform (FFT): dj = tt'L'[Dij] . Finally, the values 
resulting from Eq. 25 are smoothed before they are plotted, by 
averaging over 3 neighbouring array elements using a top-hat 
function. 

For point-masses in quasi-circular motion one obtains the 
relation (e.g. Cutler & Flanagan, 1994) 



d£P 
d/ 



G 



1/3 



2 

TT C 



with fo = -— T — 
123 p 



(26) 



The energies given in Eqs. 26 and 26 represent averages over 
one orbital revolution. 



2.4. Neutrino treatment 

In the context of the three-dimensional hydrodynamic model- 
ing of the coalescence of binary neutron stars, the production 
and emission of neutrinos is treated in terms of a neutrino- 
leakage scheme. Since the timescale of the dynamical merging 
event covered by the presented simulations is of the order of 
10 ms only, transport effects by the diffusion of neutrinos can 
be expected to be of minor importance. Moreover, neutrino mo- 
mentum transfer plays a negligible dynamical role. Therefore 
the use of a leakage scheme that takes into account the possi- 
ble equilibration of neutrinos and the opaque, hot stellar matter 
appears adequate to describe local neutrino sources and sinks. 

Electron-type neutrinos Vg and Dg are created via charged- 
current /3-processes 



+ p — > n + Vg , 
+ n — > p + Dg , 
by electron-positron pair-annihilation, 

e~ + — > Vi + Vi , 



(27) 
(28) 

(29) 



and - with a dominant rate at high densities and high electron 
degeneracy - by plasmon decays. 
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Vi + Vi 



(30) 



The emission of Vg and Vg extracts lepton number and energy, 
while heavy-lepton neutrinos = p^, v^, Vt, i't, generated 
viathe "thermal" processes of Eq. (29)andEq. (30), carry away 
energy from the stellar medium. At low optical depths the pro- 
duction and emission of neutrinos is computed directly from 
the rates of the above processes. In contrast, at high optical 
depths the equilibration timescales are much shorter than the 
timescales of neutrino diffusion and of hydrodynamic changes. 
Neutrinos are therefore assumed to be present in their chem- 
ical equilibrium abundances and the loss of neutrino number 
and energy proceeds on the diffusion timescale. For detailed 
information about the neutrino production rates and the super- 
position of the treatments in optically thick and thin regions, see 
Appendix B. All rates are approximated by simple analytical 



formulae, which allow for a very fast numerical evaluation and 
whose accuracy of the order of a few 10% in describing the 
neutrino emission of hot neutron star matter is sufficient in the 
context of the presented work. 

The optical depth for Vg and Vg is dominated by the in- 
verse reactions of the /3-processes of Eq. (27) and Eq. (28), Vg 
absorption onto neutrons. 



n + Vg 



+ P, 



and Vg absorption onto protons, 

p + Vg — > + n . 



(31) 



(32) 



An important contribution also comes from neutral-current scat- 
tering off nucleons. 



Vi + 



(33) 



which is by far the most important opacity source for heavy- 
lepton neutrinos. Since for numerical reasons the employed 
equation of state of Lattimer & Swesty (1991) had to be used 
as a table (see Sect. 2.2) which we desired to keep as small 
as possible, we did not tabulate detailed information about the 
baryonic composition of the neutron star medium. Therefore 
we neglected neutrino scatterings off nuclei and assumed the 
stellar gas to be completely dissociated into free neutrons and 
protons in the determination of optical depths. The associated 
uncertainties are of the order of a few 10% and therefore within 
the linoits of accuracy of other aspects of the neutrino description 
in this work. The spectrally averaged opacities for the processes 
of Eq. (3 1)-Eq. (33) are given in Appendix A. 

Once the neutrino opacities are determined one can com- 
pute optical depths and diffusion timescales, which then enter 
the computation of the "effective" lepton number and energy 
loss rates of the stellar gas. The latter are defined as a continuous 
superposition of the rates of the processes Eq. (27)-Eq. (30) at 
low optical depths and of a diffusion loss term that represents 
the neutrino leakage from optically thick regions on the neutrino 
diffusion timescale (see Appendix B). Adding up the neutrino 
emission from all zones of the grid yields total neutrino num- 
ber and energy fluxes. This neutrino leakage treatment was 
compared with the results of neutrino diffusion calculations for 
cooling (spherically symmetrical) protoneutron stars and the 
agreement was found to be on the level of a few 10% again. 

3. Initial conditions 

The initial situation is defined by two identical neutron stars 
of about M = 1.63 Mq and 15 km radius that are placed at a 
center-to-center distance of ao = 42 km on a grid of size 82 km. 
The distributions of density p(r) and electron fraction Yg(r) 
are taken from a one-dimensional (cold) neutron star model in 
hydrostatic equilibrium. The radius of 1 5 km is obtained for the 
Newtonian case. For a cool neutron star with baryonic mass of 
M = 1 .63 Mq the general relativistic stellar structure equations 
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Fig. 2. Density, temperature, electron fraction, and enclosed baryonic 
mass as functions of radius for the initial neutron star model. 



yield a gravitational mass of about Mg « 1 .5 Mq and a radius 
of 11.2 1cm with the equation of state of Lattimer & Swesty 
(1991). A softer equation of state than the one used in this work 
would imply a more compact star and a lower gravitational mass 
for the same baryonic mass. 

Note that we do not start our simulations with cold, T = 0, 
initial models. Due to the extreme sensitivity of the temperature 
to small variations of the internal energy in the degenerate limit 
(see Fig. 1), we decided to set the initial temperatures in our 
neutron stars to small, but finite values. By choosing temper- 
atures corresponding to thermal energy densities of the order 
of 3% of the minimum internal energy density for a given den- 
sity p and electron fraction (which is the degeneracy energy 
density at T = 0) we can avoid a great part of the tempera- 
ture fluctuations. These are associated with the inversion of the 
equation of state for the temperature and are caused by small 
numerical errors in the internal energy as it is computed as the 
difference of total and kinetic energies in the PPM scheme. 
Fig. 2 shows density, temperature and electron number fraction 
of the initial neutron star model. The two-dimensional distribu- 
tion of the initial temperature in the orbital plane can be seen in 
panels b of Figs. 4, 6, 14 and 16. 

During the run, non-physical temperature spikes in isolated 
zones appear mainly due to the low spatial resolution combined 
with the highly degenerate state of the neutron star matter. These 
spikes also lead to non-conservation of energy. We reduce sin- 
gular temperature spikes in each of the three axial directions 
by averaging over the two neighbouring zones. We do this only 
for zones in which the temperature was over 9 MeV, the den- 
sity less than 1.5 • 10'"* g/cm^, and the temperature contrast in 
both flanks of the spike are larger than 1 MeV. This procedure 
changes only the outermost surface layers of the neutron star 
and its negligible influence on the model evolution and the re- 
sults was corroborated by the higher resolved model A128 in 
comparison with model A64. 

The density of matter in the surroundings of the neutron stars 
is set to 10^ g/cm^ with an internal energy equal to the value in 
the neutron star at the same density. To avoid this low-density 
matter being accelerated and falling onto the neutron stars. 



the velocities and kinetic energy are reduced to zero in zones 
in which the density is less than 3 • 10^ g/cm^. Additionally, 
the temperatures are also reset to their initial values where the 
density is smaller than 4 • lO'" g/cm^. These changes only affect 
relatively low-density material in the computational volume not 
filled by neutron star matter. 

With a numerical grid of roughly 1 km spatial resolution one 
cannot accurately represent the density decline near the surface 
of the neutron stars, where the density typically decreases by 
a factor of ten (one dex) every 100 m over a distance of about 
700 m. We artificially soften the edge by imposing a maximal 
change of 2 dex from zone to zone. The thickness of the surface 
layer thus results to 3 zones. 



Table 1. Parameters and some computed quantities for all models. 
JV is the number of zones per dimension in the orbital plane, S is 
the direction of the spin of the neutron stars, Tmax is the maximum 
temperature reached on the grid, £ is the maximum gravitational-wave 
luminosity, £ is the total energy emitted in gravitational waves, h is 
the maximum ampUtude of the gravitational waves at r = IGpc. 
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We prescribe the orbital velocities of the coalescing neutron 
stars according to the motions of point masses, as computed 
from the quadrupole formula. The tangential components of 
the velocities of the neutron star centers correspond to Kepler 
velocities on circular orbits, while the radial velocity compo- 
nents reflect the emission of gravitational waves leading to the 
inspiral of the orbiting bodies (see Kokkotas & Schafer, 1995): 



Vt = ujr = 



Vr = r = 



GM 



2ao 
16 



(34) 
(35) 



for r = ao/2 and equal neutron star masses. Eq. 35 is valid in 
ADM coordinates (in harmonic coordinates the factor 16 is to 
be substituted by a factor 64). The initial angular frequency of 
the Newtonian Keplerian orbit is 



2GM 



U) = 



(36) 



and the period follows to P = In/ ui = 2.6 ms for ao = 42 km 
and M = 1.63Mq. Even if the neutron stars are born on very 
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elliptic orbits, gravitational-wave emission drives the eccentric- 
ity towards zero before the neutron stars come close enough for 
mass exchange (e.g. Imshennik & Popov, 1994). 



A: point mass 




B: solid body C: anti spin 

Fig. 3. A sketch of the three different initial spin distributions of mod- 
els A, B, and C. The circular arrows between the neutron stars (circles) 
indicate the orbital motions, while the smaller circular arrows above 
show the additional spins of the neutron stars about their respective 
centers. The straight arrows inside the circles sketch the velocity dis- 
tributions in the neutron stars: in model A all parts of the neutron stars 
are given the same velocity, while the additional spins in model B 
result in a solid body type rotation of the neutron stars. In model C the 
neutron stars are rotating opposite to the orbital spin direction. 

A spin of the neutron stars around their respective centers 
is added and varied from model to model. Table 1 lists the 
parameters and some computed quantities for all models. The 
A models do not have any additional spin added on top of 
their orbital velocities. In this case all parts of the neutron 
stars start out with the same absolute value of the velocity. In 
models B64 and C64 a spin is added for each neutron star. 
The angular velocity (both magnitude and direction) of the 
spin in model B64 is equal to the angular velocity of the orbit. 
This results in a soUd body type motion for model B64. In 
contrast, the direction of the angular velocity in model C64 is 
opposite to the orbit's angular velocity. The different cases are 
sketched in Fig. 3. We consider these three models to delimit 
in some sense extreme cases for the interaction of coalescing 
binary neutron stars. While the solid body case B64 represents a 
merging event with Uttle internal shearing motions, model C64 
is characterized by large velocities in opposite directions at the 
contact point of the neutron star surfaces and the coalescence 
is accompanied by violent shearing interaction. Model A128 
has the same velocity distribution as model A64. However, it 
is computed with double the resolution in all dimensions and 
thus yields a test for the accuracy of the 64^ simulations. The 
two-dimensional velocity distributions in the orbital plane for 
models A64, A128, B64 and C64 can be seen in panels panels a 
of Figs. 4, 6, 14, and 16, respectively. 

Model T64 is set up with the same initial conditions as 
model A64 except for the distance between the neutron stars 



which is chosen to be 60 km instead of 42 km (and the size of 
the grid consequently enlarged to 105 km a side). The larger 
initial separation of the neutron stars in model T64 corresponds 
to a time about 23 ms before two point masses reach the ini- 
tial center-to-center distance of the neutron stars in models 
A64, A128, B64, and C64. This model is run for approximately 
4.8 ms, which is about 1.1 orbital periods of two point masses 
at the given distance. Model T64 is intended to test whether the 
orbit is numerically stable for large initial separations. More- 
over, the finite extension of the neutron stars leads to differ- 
ences compared to the point-mass case and internal friction is 
checked in its effect on the temperature evolution during the 
spiral-in phase. 

In Eq. 9 the term +2P^^ can be replaced by the analyti- 
cally equivalent term —2vi^. This can be seen by applying 
the product rule to rewrite the derivatives and using the fact 
that surface terms are neghgibly small. A test confirms that 
this replacement also does not produce any significant numeri- 
cal difference in the resulting gravitational-wave luminosity C 
(Eq. 5). When model A64 is rerun for approximately 1 .5 ms the 
difference of C between the two cases does not exceed 0.5% 
during this time. 

The initial distance, at which the neutron stars are placed 
at the begiiming of the simulation, is a compromise between 
physical accuracy and computational resources. On the one 
hand we would like to start with the neutron stars at a large dis- 
tance in order to correctly simulate the tidal deformation they 
experience during inspiral. However, since the rate at which 
gravitational waves radiate energy away from the system in- 
creases with the fifth power of the distance (Eq. 24) the time to 
coalescence (which is proportional to the fourth power of the 
distance, Eq. 20) increases by a huge amount, too. The test cal- 
culations with an initial separation of the neutron stars of 60 km 
reveal a negligible orbital decay of about 1% radius decrease 
during 1 ms (one quarter orbit). Starting all simulations at such 
a large distance would be numerically prohibitive. 

On the other hand, tidal deformation studies have been done 
(e.g. Lai et al, 1994; Reisenegger & Goldreich, 1994) which 
yield tidal deformations for extended objects of around 20% 
(for the principal axis) and normal mode excitation of at most 
a few percent of the radius at a distance of about 2.8 radii. 

Two more test calculations clarify the influence of gravi- 
tational wave emission and of the dynamical tidal instability 
(Lai et al, 1994, and references cited therein). In one model the 
gravitational wave emission (and backreaction) is switched off 
and the initial orbit is taken to be circular. All the other ini- 
tial conditions remain as described (in particular the separation 
of 42 km). In this case the neutron stars also merge, however 
on a longer timescale of approximately 2.5 revolutions. This 
merging is due to the fact that the two neutron stars are not in 
equiUbrium in their common gravitational field at the begiiming 
of the simulations and thus start to develop tidal deformations 
and oscillations which eventually bring them to the point of 
dynamical instability. Note that the initial separation of the two 
neutron stars (2.8 km) is only slightly larger than the mini- 
mum distance (2.6 km) for dynamical instability of objects in 
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equilibrium (Lai et al, 1994). However, without gravitational 
wave emission the coalescence proceeds on a timescale that is 
2-3 times longer. We therefore consider the approximation of 
spherical neutron stars instead of an initial equilibrium config- 
uration as acceptable. In a second test model the gravitational 
wave emission (and backreaction) is again switched off, but the 
initial orbital velocities are given a radial component (Eq. 35) 
as if the preceding evolution had been influenced by gravita- 
tional wave emission. In this case the two neutron stars merge 
on nearly the same timescale as in the simulations that include 
gravitational radiation. The reason is the following. Since we 
start our calculations with the two neutron stars close to the dis- 
tance where dynamical instability sets in, any additional radial 
velocity that decreases the orbital separation quickly pushes the 
neutron stars beyond the linnit of dynamical instabiUty. 

Viscosity determines the velocity distribution within the 
neutron star: if it is large enough, spin-up during inspiral leads 
to tidal locking. However, Kochanek (1992), Bildsten & Cutler 
(1992) and Lai (1994) argue that it is very unlikely that micro- 
scopic shear and bulk viscosity alone are sufficient for spin-up 
or for tidal heating of the neutron stars to temperatures of more 
than 10^ K. For this reason, we consider model A that contains 
two non-corotating neutron stars as generic model. Models B 
and C represent variations obtained by adding spins to the in- 
dividual neutron stars. 

The calculations of models A64, B64 and C64 were per- 
formed on a Cray-YMP 4/64, needed about 16 MWords of main 
memory and took approximately 40 CPU-hours per model. 
Model A128 was computed on a Cray-EL98 4/256, required 
about 22 MWords of memory and 1700 CPU-hours. 

4. Numerical results 

4.1. Dynamic evolution 

4.1.1. Models A64, A128 and T64 

We start our description of the hydrodynamic merging events 
with the models A64 and A128, which basically differ only in 
the resolution, model A 128 having double the number of zones 
per dimension than model A64. Figs. 4 and 5 show the den- 
sity, velocity, and temperature for the temporal evolution of 
model A64, and Figs. 6 and 7 give the corresponding informa- 
tion for model A 128. The initial distance of the neutron stars 
in our models (approximately 2.8 neutron star radii) is very 
close to the separation where the configuration becomes dy- 
namically unstable, which is approximately at a distance of 2.6 
radii (Lai et al, 1994, and references cited therein). Therefore 
it is not surprising that already after one quarter of a revolu- 
tion (ftj 0.6 ms, cf. Eq. 36) the neutron star surfaces touch due 
to the tidal deformation of the stars (Panels 4c and 6c). The 
parts of the neutron stars closest to the orbital axis experience 
the strongest tidal stretching and their temperatures decrease 
below 2 MeV (Figs. 4d and 6d), while the temperature of the 
bulk of the matter towards the geometric centers of the neutron 
stars falls from the initial 6 MeV to 4MeV. 



As soon as the surfaces of the two neutron stars touch, the 
kinetic energy of the matter that moves in opposite directions 
is converted into internal energy and an elongated bar of hot 
material forms. Two fluid vortices appear at this time (Fig. 4e 
and 6e) and are associated with two separate hot "spots" (Fig. 4f 
and 6f). Isolated regions of high temperature can be observed 
during the whole subsequent evolution. Due to the strong gravity 
field of the massive, merging object the pressure varies smoothly 
within the merged object. Therefore regions with high temper- 
ature are associated with a lower density. This is visible in all 
plots when comparing the temperature and density contours, 
e.g. in Fig 6e and 6f. 

Figure 8 shows the separation between the density maxima 
(not centers of mass) of the neutron stars. One notices, that for 
model A64 the neutron stars merge within 1 .5 ms and then two 
density maxima form and separate again after about 2.5 ms. 
This swinging of the neutron stars, caused by their angular mo- 
mentum, is not visible for model A128 in this figure. However, 
comparing the contours of the density in Figs. 5a and 7a, one 
recognizes that also in model A 128 the two initial neutron stars 
can still be discriminated. 

The evolution of the temperature distribution differs be- 
tween models A64 and A 128 after 4 ms as can be seen from a 
comparison of Figs. 5d and f with Figs. 7d and f, and is also 
clear from Fig. 9. This must partly be caused by the dissolving 
of the two fluid vortices, which are associated with the spots 
of high temperatures. In Fig. 5c the vortices are still present, 
while they are nearly completely dissolved and broken up into 
smaller ones in the better resolved model A 128 at around the 
same time. Fig. 7c. The maximum temperature on the grid is 
shown in Fig. 9. The temperature increases by roughly 20 MeV 
for the A models when the surfaces of the two neutron stars 
come into contact. 

However, because the matter is degenerate and the pressure 
not very sensitive to the temperature, the evolution of the den- 
sity distribution is fairly similar in model A64andmodel A128: 
Panels e in Figs. 5 and 7 show one merged object that contains 
most of the mass and is surrounded by a distended thick disk. 
In both models A64 and A128 the density contour correspond- 
ing to p = lO'^ g/cm-' is approximately circular and extends 
nearly to the edges of the grid, i.e. has a diameter of approxi- 
mately 80 km. The maximum density on the grid is displayed 
in Fig. 10 and is marginally higher for model A 128 compared 
to the models with a resolution of 64 zones per dimension 
(7.2 ■ lOi'^ g/cm3 compared to 6.9 ■ lO^'' g/cm^). 

The merged object has a total mass of approximately 3 Mq, 
which is larger than the maximum mass of stable neutron stars 
for the equation of state of Lattimer & S westy ( 1 99 1 ) in the gen- 
eral relativistic case. Therefore it is of interest to compare the 
size of the object with its Schwarzschild or event horizon. Since 
our computations are basically Newtonian it is not easy to es- 
timate possible implications of a general relativistic treatment. 
A problem arises from the difficulty to choose the appropriate 
coordinate system for the discussion. We shall refer to har- 
monic coordinates or ADM (generalized isotropic) coordinates 
which are adequate for binary systems. Coordinate radii f in 
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Fig. 4. Contour plots of model A64 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for the 
temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dcx, while the temperature contours are linearly 
spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm'. The legend 
in each panel at the top right comer gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 
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Fig. 5. Contour plots of model A64 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for the 
temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dcx, while the temperature contours are linearly 
spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm'. The legend 
in each panel at the top right comer gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 
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Fig. 6. Contour plots of model A128 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for 
the temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dex, while the temperature contours are 
linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm^. 
The legend in each panel at the top right corner gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 
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Fig. 7. Contour plots of model A128 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for 
the temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dex, while the temperature contours are 
linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm^. 
The legend in each panel at the top right corner gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 
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Fig. 8. The separation of the density maxima of the two neutron stars as a function of time for all five models. The thin, monotonously 
decreasing curve shows the distance between two inspiraling point-masses (Eq. 20). Note, that two point-masses at an initial distance equal to 
the center-to-center distance of the neutron stars of model T64 need about 23 ms to reach the initial separation of models A64, A128, B64, 
and C64. 
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Fig. 9. The maximum temperature on the grid as a function of time for all five models. 



Schwarzschild coordinates, f in harmonic coordinates, and r* 
in isotropic coordinates are related by the following equation: 

r = r + 0.5Rs = r*il + 0.25R,/r*f , (37) 

where = 2GM/c^ is the Schwarzschildradius of a field 
generating mass M. The function m(r) which gives the mass 
of the merged object enclosed by radius r is shown in Fig. 1 1 . 
For a given mass, the radius rim) at which the object attains 
this mass, is always larger than 0.5 i?s(m(r)), the radius of the 
Schwarzschild horizon in harmonic coordinates. (i?s(m(r)) = 
2Gm{r)l(? is the Schwarzschildradius that corresponds to the 
enclosed mass m(r).) However, for the mass range between 
about 1 and 3 Mq the radius 3 Rg lies outside the merged object. 



3Rs is, with an uncertainty of roughly 25%, suggested to be the 
radius of the last stable circular orbit for an equal-mass binary 
in harmonic coordinates (Kidder et al, 1992; Wex & Schafer, 
1993). This, of course, means that deeper investigations are 
needed. In Fig. 1 1 the straight hne denoted by 2.5 R^ gives 
the mass-radius relation in harmonic coordinates for the last 
stable circular orbit of a test particle with nonzero rest mass in 
a Schwarzschild field. In ADM coordinates the radius of the 
Schwarzschild horizon takes the value 0.25 Rg and the radius 
of the last stable circular orbit of a nonzero-mass test body is 
about 2.475 Rs- The radius of 3.0 R^ in harmonic coordinates 
is located at approximately 2.979 R^ in ADM coordinates. 
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Fig. 10. The maximum density on the grid as a function of time for all five models. 
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Fig. 11. The bold solid, dashed and dotted curves show the enclosed 
mass m(r) as a function of radius (starting at the center of the coalesced 
object, i.e., at the grid center) at the end of the simulations. The three 
straight lines depict multiples of 0.5, 2.5, and 3.0 of the Schwarzschil- 
dradius Rs(m{r)) that corresponds to the mass on the vertical axis. In 
harmonic coordinates, 0.5 Rs defines the radius of the Schwarzschild 
(event) horizon, 2.5 Rs is the smallest radius of stable circular orbits 
for test particles with nonzero rest mass in a Schwarzschild field, and 
3 Rs denotes (roughly, not better than about 25%) the closest radius 
of stable circular orbits for an equal-mass binary. 



When the merged object rotates so fast that centrifugal 

forces are important the collapse to a black hole might be de- 
layed by the stabilizing effect of these forces. Fig. 12 compares 
the Kepler velocity at a given radius, which is only dependent 
on the mass inside that radius, to the average model velocities at 
the same radius. One can see that the actual velocities are some 



Fig. 12. The upper curves show the Kepler velocity Vk{m(r)) that 
corresponds to the mass inside radius r as seen in Fig. 11. The lower 

curves depict the average orbital velocity in the merged objects as 
function of the distance from the grid center. 



factor of 4 smaller than the Kepler velocities, so centrifugal 
forces are not expected to be strong enough to prevent or delay 
the collapse to a black hole. 

The density and temperature distributions of the test 
model T64 can be found in Fig. 13 for two snapshots in time. 
The first two panels a and b show the initial state, while the 
panels c and d show the distributions after about 1 .2 orbital 
revolutions. Only relatively minor changes of the structure oc- 
cur during this time: the thickness of the surface layers varies, 
the shapes of the neutron stars are slighly deformed from the 
original spheres, and the internal temperatures in the stars de- 
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Fig. 13. Contour plots of model T64 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for 
the temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dex, while the temperature contours are 
linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm?. 
The legend in each panel at the top right comer gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 



crease due to the net expansion and increase of the volumes. The 
maximum density (shown in Fig. 10) oscillates with a period 
of approximately 0.3 1 ms (which is the sound crossing time of 
the neutron stars) around an average value of 5.6 • 10'"* g/cm^. 
This oscillation has a damping timescale of approximately 
2 ms. Note that the maximum density of the initial neutron star 
model is 5.7 ■ lO''' g/cm^. The maximum temperature fluctu- 
ates around 8 MeV (Fig. 9), but these temperatures affect only 
few zones none of which appear in Fig. 13. These results for 
test model T64 show us that it is rather unhkely that the evo- 
lution preceding the initial configurations of our models A64, 



A128, B64, and C64 will produce strongly distorted and heated 
neutron stars. 

4.1.2. Models B64andC64 

Model B64 (Figs. 14 and 15) represents the evolution of two 
neutron stars that rotate like one rigidly rotating solid body. 
The corresponding initial velocity distribution was adopted in 
a number of other investigations, e.g. in most of the models 
computed by Shibata et al (1993) and Nakamura & Oohara 
(1991), and in earlier works of this group. The assumption of 
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Fig. 14. Contour plots of model B64 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for 
the temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dcx, while the temperature contours are 
linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm^. 
The legend in each panel at the top right comer gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 
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Fig. 15. Contour plots of model B64 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for 
the temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dcx, while the temperature contours are 
linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm^. 
The legend in each panel at the top right comer gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 
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Fig. 16. Contour plots of model C64 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for 
the temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dex, while the temperature contours are 
linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm^. 
The legend in each panel at the top right comer gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 



20 



M. Ruffert et al.: Coalescing neutron stars - a step towards physical models 



Fig. 17. Contour plots of model C64 showing cuts in the orbital plane for the density together with the velocity field (left panels) and for 
the temperature (right panels). The density contours are logarithmically spaced with intervals of 0.5 dcx, while the temperature contours are 
linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values, the density is measured in units of g/cm^. 
The legend in each panel at the top right comer gives the scale of the velocity vectors and the time elapsed since the beginning of the simulation. 
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rigid rotation appears natural when tidal locking is assumed due 
to the action of viscous forces. This situation, however, does not 
appear to be likely (see discussion in Sect. 3). 

Model B64 shows important differences in the initial veloc- 
ity distribution when compared with model A64. The velocities 
at the neutron star surface regions that face each other are much 
smaller in model B64. When these surfaces touch, the dissipa- 
tion of kinetic energy into heat is much weaker. This can be 
seen very clearly by comparing Fig. 14f with Fig. 4f. While in 
model A64 the temperatures have reached 14 MeV after about 
1.7 ms, they are only 8 MeV in model B64. This difference 
is also visible in the maximum temperatures (Fig. 9), where 
model B64 is the only model (apart from the test model T64) 
that does not show a sudden temperature jump at approximately 
1 ms. Although eventually also in model B64 two hot spots de- 
velop (Fig. 15d), these are less prominent and dissolve towards 
the end of the simulation (Fig. 15f) and form a ring of hot 
material in the orbital plane. Again (cf. model A128), this dis- 
solution happens at the same time when the vortex motions 
disappear. 

Differences of the velocity distribution of models A64 
and B64 also exist at those parts of the surfaces of the neutron 
stars that are most distant from the orbital rotation axis. Since 
the velocities of model B64 are much larger here, the centrifu- 
gal forces are also much stronger, which in turn facilitates mass 
shedding from the outer layers. Thus, two prominent spiral arms 
form (Fig. 15a) from which matter is lost off the grid. Fig. 18 
shows the amount of gas that flows off the grid as a function of 
time for all models. Model B64 sheds approximately five times 
more matter than model A64, and model A128 is another factor 
of two below model A64. Although initially prominent, the spi- 
ral arms (on our grid) do not seem to survive several revolutions 
and after about 4.5 ms a thick disk, just as in model A64, has 
developed. Fig. 15e does not differ quaUtatively from Fig 5e or 
Fig 7e. 

Most of the mass that flows off the grid remains bound 
to the system. Fig. 19 shows the cumulative amount of matter 
that gets unbound as a function of time. For every zone along 
the grid boundaries, we check whether the sum of potential, 
kinetic, and thermal energy of the outward streaming gas is 
positive. To estimate the amount of unbound matter we add up 
the mass flow across the boundaries only for these zones. We 
find that a few times 1O~"*M0 might be able to escape from 
the system. These numbers are obtained by using the thermal 
energy of the gas which is the differerence between its internal 
energy and its internal energy at T = for given values of 
density and Yg . Taking the total internal energy instead of just 
the thermal energy could lead to somewhat larger estimates of 
the mass that gets gravitationally unbound. In case of model B64 
the difference is about a factor of 3 (Fig. 19). Neutrino-energy 
deposition in the cool outer regions of the disk could power 
a low-density mass flow like the neutrino-driven wind from a 
newly formed neutron star in a Type II supernova. Moreover, 
energy release from the recombination of nucleons or from 
nuclear reactions (Davies et al, 1994) could have dynamical 
effects and could increase the mass loss from the disk. All these 



processes suggest that our numbers for the unbound mass are 
only lower estimates. The stronger gravitational potential in a 
fully general relativistic treatment, on the other hand, would 
imply a stronger binding of the disk material and would lead to 
a corresponding reduction of the mass loss. 
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Fig. 18. The cumulative amount of matter that flows off the grid during 
the simulations as a function of time for models A64, A128, B64, 
and C64. 



If one interprets model B64 as emerging from model A64 
by adding spins to the neutron stars around their respective axes, 
then model C64 is constructed from model A64 by adding neg- 
ative spins. Thus, contrary to model B64, the velocities in C64 
are largest at points closest to the orbital rotation axis (Fig. 16a). 
The differences to model A64 as discussed for model B64 
above are now found to be present in the reverse direction. As 
soon as the neutron star surfaces touch each other the temper- 
atures quickly rise to high values, roughly 30 MeV (Fig. 16f), 
and later to even more than 45 MeV (Fig. 9). As already seen 
in model A64, the hot contact region of the neutron stars splits 
into two isolated hot spots which remain present until the end 
of the simulation (Fig. 17b, d, and f). 

The high velocities at the neutron star surfaces facing each 
other also lead to a stronger deformation of the approach- 
ing stars in model C64 (Fig. 16c) compared to models A64 
(Fig. 4c) and B64 (Fig. 14c), the latter two developing ellip- 
soidal deformations with the major axes of the neutron stars be- 
ing aligned shortly before coalescence. On the other hand, the 
velocities of the parts of the neutron star surfaces facing away 
from the orbital rotation axis are small and no pronounced spi- 
ral arm extensions like in model B64 appear (compare Fig. 16e 
with Fig. 14e and Fig. 17a with Fig. 15a). 

At the end of our simulations (Figs. 5e,f, 15e,f, 17e,f) the 
coalesced object has the following structure. Most of the mass 
(96%) is contained in a high-density central, rotating body that 
is slightly triaxial. The major axes of the isodensity contours 
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Fig. 20. Contour plots of the four models A64, A128, B64, and C64 showing the density distribution in cut planes perpendicular to the orbital 
plane. The density contours are logarithmically spaced with intervals of 0.5 dex. The bold contours are labeled with their respective values. The 
times at which these snapshots are taken are the same as the times of panels f of Figs. 5, 7, 15, and 17, respectively. 
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Fig. 21. Contour plots of the four models A64, A128, B64, and C64 showing the temperature distribution in cut planes perpendicular to the 
orbital plane. The temperature contours are linearly spaced with 2 MeV increments. The bold contours are labeled with their respective values. 
The times at which these snapshots are taken are the same as the times of panels f of Figs. 5, 7, 15, and 17, respectively. 
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Fig. 19. The cumulative amount of matter that is unbound when flowing 
off the grid as a function of time for models A64, B64 and C64. The 
thick dashed, dotted, and solid lines, respectively, are computed with 
the criterion that the sum of gravitational, kinetic, and thermal energy 
of the gas are positive. In this case no matter gets unbound in the better 
resolved model A128. Using the total internal energy instead of only 
the thermal energy leads to the result given by the thin solid line for 
model B64 (denoted by B64*). 



are not aligned. The central body wobbles and rings and contin- 
uously sends pressure and density waves into the surrounding, 
thick, distended disk which is formed by the matter in the outer 
parts of our computational grid. These waves transfer angular 
momentum and influence the disk's structure and extension, 
in particular cause periodic expansion and contraction motions 
with spiral-wave-like appearance. Moreover, the interactions 
of the waves generates heating of the disk gas over the sim- 
ulated period of several milliseconds after the merging of the 
neutron stars. Probably also most of the material that flows off 
the grid (about 0.1 Mq) would add to the disk, because we 
expect only very little (a few 10~'*Mq) of this matter to get 
unbound. Maximal temperatures occur in hot spots located in 
a ring of material at a distance of about 10-15 km around the 
central density maximum. The densities in this very hot layer 
are typically l-2-10''*g/cm^. 

Density and temperature distributions for models A64, 
A 128, B64, and C64 in two orthogonal planes perpendicular 
to the orbital plane are shown in Figs. 20 and 2 1 . The displayed 
snapshots correspond to the times of the last panels f of Figs. 5, 
7, 15, and 17, respectively. The high-density cores are essen- 
tially spherically shaped in all of the cases of Fig. 20, only in 
model A128 small cuspy bulges of the p = lO"' g/cm^ contour 
at the equator (orbital plane) indicate the very fast rotation of the 
massive central object. Also, model A 128 has steeper density 
gradients at the poles than model A64, which suggests that the 
vertical extension of the disk and its shape are somewhat depen- 
dent on the numerical resolution. With higher numerical resolu- 



tion one should expect a more compact and more elliptically de- 
formed disk structure vertical to the orbital plane. Asymmetries 
of the disk to both sides of the center and differences between 
the two cut directions result from the perturbations and waves 
sent into the surrounding gas by the dynamical contractions and 
expansions of the central object. The temperature distributions 
of Fig. 21 reveal a rather complex structure. The ring-like region 
of highest temperatures (T ^ 15-20 Me V) around the cooler 
central core (see above) has a crescent or banana-shaped cross 
section. Hot knots outside the equatorial plane can occur and 
are ordered in a quadrupole type manner. In all models the re- 
gions at the poles remain comparatively cool with temperatures 
similar to those present in the inner core (T ^ 15-20 Me V). 
The outer parts of the disk have temperatures around 2-4 MeV 
in regions where the density is a few lO'^ g/cw? or less. 

Since our calculations are Newtonian, the merged object 
with a mass of about 3 Mq remains stable and does not collapse 
into a black hole. However, in a general relativistic treatment 
gravitational instability would be unavoidable on a timescale of 
a few milliseconds from the start of our modelling, because the 
mass of 3 Mq is beyond the maximum stable neutron star mass 
(about 2.2 Mq) associated with the used equation of state. 
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Fig. 22. The total emission rate of gravitational wave energy as a func- 
tion of time for models A64, A128, B64, and C64. The monotonously 
rising curve shows the emission for a point-mass binary (Eq. 24). 



4.2. Gravitational waves 

The gravitational wave luminosity (Eq. 5) is shown as a function 
of time in Fig. 22. Initially, the luminosity of our models is 
nearly equal to the emission of binary point-masses. After about 
1 ms the dynamical instability sets in and the neutron stars 
merge. Therefore, within less than 1.5 ms the power of the 
models reaches a very high, narrow maximum. Model A has a 
second, less pronounced luminosity peak at about 3.5 ms, while 
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Fig. 23. The integrated energy emitted in gravitational waves as a 
function of time for models A64, A128, B64, and C64. 



model B shows a second, broad hump around 5 ms. Model C 
is exceptional because its first luminosity maximum is only 
roughly half as high as that of the other models and is therefore 
only slightly stronger than the following, second maximum at 
about 2.75 ms. 

There seems to be a trend from model B over A to C, as- 
sociated with the decrease of the total angular momentum (see 
Fig. 3): the relative height of the first to the second luminosity 
maximum decreases. 

The time integral of the luminosities, i.e. the energy emitted 
in gravitational waves, is shown in Fig. 23. Although the main 
phases of emission occur between 1 and 4 ms, a slow rise of the 
curves towards the end of the simulations indicates that even 
during the late stages of the simulated evolution a considerable 
amount of energy continues to be emitted. The total energy 
radiated by our models during the computed neutron star coa- 
lescence and subsequent evolution lies between 1 .5% and 2.2% 
of Mqc^. 

In Figs. 22 and 23 differences of the gravitational wave 
emission between models A64 and A128 are clearly visible. 
Model A64 has a much stronger second luminosity peak and, 
correspondingly, its total energy emitted in gravitational waves 
is roughly 40% higher. These discrepancies are most likely a 
consequence of the coarser numerical resolution of model A64. 
Although the structure of the PPM code explicitly conserves en- 
ergy and momentum (in the absence of a gravitational potential 
term), this does not hold true for the angular momentum. During 
the early, very dynamical phase of the merging of the neutron 
stars, fluid whirls and vortices form which contain a significant 
fraction of the kinetic energy and angular momentum of the 
initial orbital motion. The distribution and amount of angular 
momentum plays an essential role for the subsequent behavior 
and structure of the merged object. This is obvious from the 
comparison of models A64, B64, and C64. During the merging 



process minor deviations between models A64 and A 128 are 
already visible at t « 0.7 ms (Figs. 4c,d and 6c,d), increasing 
differences at i « 1 .7 ms (Figs. 4e,f and 6e,f), and significant 
ones at t X 3.1-3.3 ms (Figs. 5a,b and 7a,b). In model A128 
structures and density and temperature gradients are not only 
sharper and better represented, but, most important, the fluid 
motion appears less coherent and the flow pattern decays into 
a larger number of smaller vortices during the following evolu- 
tion. 

These differences of the fluid flow are reflected in the time 
evolution of the separation of local density maxima as displayed 
in Fig. 8. The coherence of the flow in model A64 is associated 
with a sizable amount of angular momentum in the large-scale 
flow and leads to a transient re-formation of two distinct den- 
sity maxima after about 2.5 ms and after an intermediate pe- 
riod where the neutron stars had already merged into a single, 
very compact body. In contrast, in model A 128 much angular 
momentum is located in the smaller and more fine-stuctured 
vortices and the angular momentum in the overall rotation is 
not sufficient to support a partial re- separation of the merged 
neutron stars after coalescence. This is apparent in Figs. 5a and 
7a where the high-density, central body of model A64 shows 
a very elongated shape while in model A 128 it is more com- 
pact and spherical. The discussed dynamical consequences of 
the numerical resolution of both models are also responsible 
for the differences between A64 and A128 seen in Figs. 9 and 
10. Moreover, there is a positive correlation between the grav- 
itational wave luminosity (Fig. 22) and the separation of the 
density maxima in the merging system (Fig. 8): the gravita- 
tional wave emission of models A64 and A128 is very sinoilar 
up to about 2.5 ms which is the time when the two density max- 
ima of the neutron stars in model A64 start to reappear. Between 
2.5 ms and about 4.5 ms A64 shows its pronounced second lu- 
minosity peak, which perfectly corresponds to the duration of 
the intermediate phase of the re-formation of two distinct, lo- 
cal density maxima. After 4.5 ms both models A64 and A128 
have developed a nearly spherical high-density core and their 
gravitational wave emission becomes rather weak. 

A detailed evaluation of our models reveals that during the 
first 1-2 ms, which is the phase of the violent, dynamical merg- 
ing (see Fig. 8), model A64 loses about 10% of its initial angular 
momentum, while in the case of A128 the violation of angular 
momentum conservation seems to be negligible before about 
4.5 ms. Only during the subsequent evolution (t ^ 4.5 ms) 
the cumulative loss of angular momentum grows roughly lin- 
early with time in A 128 to reach a few per cent at the end of 
our simulation. This decrease of the angular momentum in the 
computational volume is primarily caused by the flow of matter 
off the grid, see Fig. 18, which occurs at a constant rate dur- 
ing the late evolution of model A 128. Model A 128 seems to 
be well resolved and angular momentum conserved to a high 
degree of accuracy. Since gravitational waves do not only carry 
away energy but also angular momentum it is necessary to dis- 
tinguish carefully between the physical losses mentioned above 
and the numerical destruction of angular momentum. To check 
our analysis we compare with a simulation of the dynamical 
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coalescence where the effects of gravitational waves are com- 
pletely switched off. In case of a computation with 64^ grid 
points the angular momentum is satisfactorily conserved until 
about 0.7 ms but between 1 and 2 ms about 6% of the angular 
momentum are destroyed. Although the merging proceeds in 
a different way when gravitational waves are disregarded, the 
similarity of the numbers suggests that our conclusions about 
the numerical loss of angular momentum should be correct. 

The gravitational waveforms that correspond to these lumi- 
nosities are shown in Figs. 24 and 25. In both figures the results 
of the numerical models are plotted by bold lines, while the an- 
alytic results for the point-mass binary are plotted by thin lines. 
For model T64 the time t « —23 ms coincides with the begin- 
ning of the simulations. The finite extension and deformation of 
the neutron stars in the test model T64 produces noticable de- 
viations from the point-mass solution after approximately one 
revolution (Fig. 25). A part of this effect, however, may also be 
caused by the low resolution of this model with only 9 mesh 
zones on the length of one neutron star radius. 

The amplitude and frequency of the waves in the point- 
mass approximation increases monotonously until the two point 
masses finally meet at t « 7 ms (Fig. 24). For an initial phase 
of about 1 ms the waveforms of the numerical models A, B, 
and C coincide well with those of the point-mass binary. As 
soon as the neutron stars get tidally deformed, the waveforms 
deviate in amplitude and frequency. The merged neutron stars 
tend to produce waves at a dominant frequency between 1.5 
and 2 KHz as can be seen in Fig. 26 which gives the Fourier 
transform of the gravitational waveforms of Fig. 24. Beyond 
the broad maximum the power spectra of the gravitational wave 
emission show a steep decline towards higher frequencies with 
an average power-law index of about —7....— 10. 

The cumulative emission of gravitational-wave energy as a 
function of frequency is shown in Fig. 27. In the upper part of 
each panel, the downward sloping straight line represents the 
energy loss per unit frequency interval of a point-mass binary 
as given by Eq. 26. The frequencies to the left of the vertical 
line correspond to the wave frequencies that are emitted before 
the start of the numerical modelUng, i.e. for times t < (see 
Fig. 24). We produce a combined wave by using the quadrupole 
moments for a point-mass binary for i < and taking the 
numerically obtained quadrupole moments (Eq. 19) at i > 
(see Zhuge et al, 1994). The combined wave is then Fourier 
analysed up to the times given in the lower left corners of 
the panels of Fig. 27. With the Fourier transform, the energy 
spectrum is calculated via Eq. 25. 

At low frequencies the energy spectrum calculated for the 
combined waves fits very well to what is expected from the 
point-mass approximation. Although the merged object in the 
simulations radiates gravitational waves for a longer time than 
the point-masses, the very much smaller ampUtudes (Fig. 24) 
result in less energy emitted at almost all higher frequencies. 
A prominent peak is visible for all models at about 2 KHz. 
By Fourier transforming the signal of the combined wave up 
to different times i > 0, we are able to roughly locate the 
time at which the peaks are produced. At around 1 .6 ms all 



Fig. 24. The gravitational waveforms, ft+(0,0) and hx(0,0}, for the 
models A64, B64, and C64 are plotted by bold, solid lines beginning at 
t = 0. The thin, solid lines represent the analytic result for a point-mass 
binary in the quadrupole approximation (Eqs. 21 and 22). 
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Fig. 27. Snapshots of the energy spectrum of emitted gravitational waves at four different times for the models A64, B64, and C64. The time 
up to which the Fourier transform is performed is indicated in the lower left comer of each panel. The straight, downward sloping line is the 
spectrum of a point-mass binary, as given by Eq. 26. The vertical lines indicate the frequency corresponding to the orbital frequency of two 
point masses at the initial distance of the neutron stars in our numerical models. 



models have a featureless, monotonously declining spectrum 
at high frequencies. Model A64 starts to produce its 2 KHz 
maximum at around 3 ms, while model B64 shows indications 
of its maximum after about 5 ms. Model C64 has developed a 
prominent 2 KHz peak already at 3 ms. 

The gravitational waves yield information about the acceler- 
ation of the tensor components of the mass quadrupole moment 
of our system or, according to Eq. 19, about the trace-free part 
of the sum of two times the system's kinetic energy tensor and 
the stress tensor for the Newtonian gravitational field. As is 
well known, taking the trace-free part is equivalent to subtract- 
ing the spherically symmetric part. Thus, the gravitational wave 
field gives only information about non-spherically symmetric 
(stress-energy) aspects of its source. 

While gravitational waves picture the (hydro-)dynamics and 
mass motions and are produced in the KHz regime, the radiation 
of neutrinos, which are emitted with thermal energies typical 



of the region around the neutrino sphere, reflects the thermo- 
dynamical properties of the emitting matter, in particular of the 
disk material that surrounds the massive central object formed 
after the merging of the binary neutron stars. The characteris- 
tics of the neutrino emission and the corresponding implications 
for gannma-ray burst models will be addressed in a forthcoming 
paper. 

5. Summary and discussion 

The hydrodynamic simulations presented in this work follow 
the density and temperature evolution as well as the gravita- 
tional wave emission of two merging neutron stars. Initially, 
each of the two identical, cool neutron stars is in hydrostatic 
equilibrium and has a baryonic mass of « 1 .6 Mq and a radius 
of « 15 km. The neutron stars are placed at an initial center- 
to-center distance of 42 km. No attempt is made to construct 
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Fig. 25. The gravitational waveforms, h+(0, 0) and hx(0, 0), for the test 

model T64 arc plotted by bold, solid lines beginning al t = —23 ms. 
The thin, solid lines represent the analytic result for a point-mass binary 
in the quadrupole approximation (Eqs. 21 and 22). 
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Fig. 26. Power spectrum (Fourier transform) of the gravitational wave- 
forms shown in Fig. 24 for the models A64, B64, and C64. 



initial equilibrium configurations of the binary system. Tlie ini- 
tial temperatures in the neutron stars are increased above the 
cold, T = 0, situation in such a way that the thermal energy 
is about 3% of the degeneracy energy for given density and 
electron fraction. This yields a central temperature of about 
8 MeV at the beginning of our simulations. The orbital veloci- 
ties of the coalescing neutron stars are prescribed according to 
the motions of point masses as computed from the quadrupole 
formula. Spins are added to the neutron stars to take into ac- 
count rotations of the neutron stars around their axes vertical to 



the orbital plane. The spins are different from model to model. 
We simulate the three cases without neutron star spins and with 
spin vectors parallel and antiparallel to the vector of the orbital 
angular momentum. The case with neutron star spins parallel to 
the orbital spin yields a rigid solid-body rotation. Recent work 
(Kochanek, 1992; Lai, 1994; Reisenegger & Goldreich, 1994) 
suggests such a small viscosity of neutron star matter that the 
stars cannot develop synchronous rotation during inspiral. This 
gives us the freedom to set up initial spins without constraints 
from the orbital motion. Contrary to Rasio & Shapiro (1994) 
who considered the corotating case only, we therefore do not 
choose our initial conditions as exact equiUbrium configura- 
tions. 

The orbit decays due to gravitational wave emission, and af- 
ter less than one revolution the stars are so close that dynamical 
instability (Lai et al, 1994, and references cited therein) sets in. 
Within 1 ms they merge into a rapidly spinning (Pspin « 1 ms), 
high-density body (p ^ 10'"' g/cm^) with a surrounding thick 
disk of material with densities p « lO'" — lO'^ g/cnP and 
rotational velocities of 0.3-0.5 c. In our Newtonian models a 
merged object with a mass of approximately 3 Mq is formed. 
This mass is far above the limiting mass of relativistic neutron 
stars for the used equation of state of Lattimer & Swesty (1991) 
and is also in excess of the Umiting masses of most currently 
discussed supranuclear equations of state. We therefore expect 
that the central object in our simulations would collapse into a 
black hole on a time scale not much longer than the dynamical 
timescale of approximately 1 ms, if the computations were per- 
formed with a fully relativistic treatment. The duration of the 
stable phase prior to collapse will depend on the masses of the 
neutron stars and on the properties of the equation of state. The 
stability of the coalesced object in the context of uncertainties 
of the supranuclear equation of state was discussed in more de- 
tail by Davies et al (1994) and the possible effects of additional 
support of the merged object by centrifugal forces associated 
with the large angular momentum was addressed by Rasio & 
Shapiro (1992) and Davies et al (1994). 

Because of the different angular momentum in the different 
models, the amount of matter that flows off the computational 
grid varies between 0.02 and 0.15 M©. Since this matter has a 
high specific angular momentum it will most likely add to the 
disk around the merged object and might even remain in the 
disk when the massive central object has collapsed into a black 
hole. At the moment when crossing the grid boundaries, only a 
fairly small fraction of this matter has the energy to escape from 
the gravitational potential of the merged binary. Our estimates 
yield a few lO"'* Mq of unbound material. Yet, momentum and 
energy transfer by pressure waves, neutrino heating, and nuclear 
energy release might lead to additional mass being stripped off 
the outer regions of the disk and could increase the mass loss. 
The use of a grid-based code like PPM does not allow us to fol- 
low the evolution of the tidally-formed spiral arms which reach 
far out from the merged neutron stars and are very prominent 
in the SPH simulations of Rasio & Shapiro (1994), Davies et al 
(1994), and Zhuge et al (1994). The formation of these spiral 
arms, however, is visible in our model with solid-body like ro- 
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tation. This model accordingly also shows the maximal amount 
of mass-shedding off the grid (0.15 Mq). The presence of spiral 
arms clearly depends on the neutron star spins and the corre- 
sponding angular momentum in addition to the orbital angular 
momentum. This conclusion can also be drawn on grounds of 
the results and figures in Shibata et al (1992). The dilute disk 
extends to the grid boundaries not only in the orbital plane but 
reaches out to the boundaries also in the vertical direction. A 
steep density gradient towards the edges of our grid and com- 
paratively high densities (lO'^-lO" g/cm^) indicate that the 
thickness of the disk certainly exceeds 40 km which is the verti- 
cal height of our computational volume. Davies et al (1994) find 
a butterfly-shaped cross section of the disk perpendicular to the 
plane of the neutron star trajectories in their SPH simulations. 

At the end of our simulations when the merger has reached 
a state near rotational equilibrium, the central, compact object 
formed after the coalescence of the neutron stars is essentially 
spherical for the model with the lowest total angular momentum 
(model C64). In case of the other models slight triaxiality occurs 
with major axis ratios of the density contour for p = lO''* g/cm^ 
of about 1.15 : 1 in the orbital plane and between 1.09 : 1 and 
1.25 : 1 perpendicular to the equatorial plane. This is in agree- 
ment with the results of Rasio & Shapiro (1994) and Zhuge et al 
(1994) who found this structure for models with stiff equations 
of state and adiabatic exponents T ^ 2.3. 

The peak luminosity of gravitational waves occurs shortly 
after the dynanoical instability of the coalescing binary has set 
in and the neutron stars start do merge. A maximum luminosity 
in excess of 10^^ erg/s is reached for about 1 ms. Although the 
energy enoission rate drops after this short but powerful out- 
burst because the merged object has a much smaller asymmetry 
than the pre-merging configuration, a significant contribution 
to the gravitational wave signal (about 50% of the total emitted 
energy), in particular at high frequencies ^ 1 KHz, is radiated 
at later times. The total energy emitted in gravitational waves 
varies between 1.7%-2.2% of Mqc^, which is considerably 
less than the values of 7-9% of MqC^ reported by Nakamura & 
Oohara (1991). The energy release in gravitational waves de- 
pends strongly on the size and compactness of the neutron stars 
since the gravitational wave luminosity increases like with 
the minimum separation of the neutron stars before they merge 
and before the binary system loses its highly non-spherical ge- 
ometry. With a radius of about 9 km the polytropic neutron star 
models of Nakamura & Oohara (1991) are much smaller than 
those of our work (about 15 km) using the equation of state of 
Lattimer & Swesty (1991). Correspondingly, the gravitational 
wave amplitudes found in our simulations reach 3 ■ 10^^^ for 
a distance of 1 Gpc, while Shibata et al (1992), like Nakamura 
& Oohara (1991) having more compact neutron stars and using 
polytropic equations of state, obtained maximum wave ampli- 
tudes of 3.5^ ■ 10^^^. Most of the gravitational wave energy is 
radiated in a window of frequencies between 1 KHz and 2 KHz 
which corresponds to the dynamical frequencies of the binary 
system and merged object. 

The phase of the maximum emission of gravitational waves 
from the time when the neutron star surfaces begin to touch 



to the moment when their density maxima merge, is accompa- 
nied by a steep increase of the gravitational wave amplitude, 
followed by an abrupt drop and a subsequent extended pe- 
riod where the radiated gravitational waves exhibit a periodic 
modulation and oscillatory behavior with indications of the su- 
perposition of modes of different frequencies. This behavior 
depends on the triaxiality of the compact, merged object and 
is caused by the wobbling and ringing of the central, massive 
body after the dynamical coalescence. The details of the mass 
flow and therefore the structure of the gravitational waveforms 
are sensitive to the amount of angular momentum in the sys- 
tem and the corresponding degree of triaxial deformation and 
asymmetry. Positive correlations between the stiffness of the 
neutron star equation of state and the nonaxisymmetry of the 
final, merged object on the one side and the strength and ampU- 
tude of the post-merging oscillations of the gravitational wave 
emission on the other were pointed out by Rasio & Shapiro 
(1994) and Zhuge et al (1994). From our simulations, aU per- 
formed with the same equation of state, we conclude that also 
the total angular momentum in the system affects the structure 
of the gravitational waveforms and determines the modulation 
of their amplitudes. In particular, the size of the gravitational 
wave amplitude, the peak luminosity of the gravitational waves, 
the strength of the post-merging emission, and the form of the 
gravitational wave energy spectrum jj, i.e., the energy emitted 
in gravitational waves at different frequencies, clearly depend 
on the rotational state of the merging system and of the coa- 
lesced object. Clear trends with the total angular momentum 
of the system are visible for some quantities, e.g., an increase 
of the maximum gravitational wave amplitude, an increase of 
the energy emitted in the 0.8-1.5 KHz band, and a decrease of 
the narrow spectral peak at a frequency of about 2 KHz. The 
temporal modulation of the waveforms and the power spectrum 
of the wave amplitude, however, seem to vary in a complex and 
non-trivial way and reflect details of the mass motions in the 
merger. 

Zhuge et al (1994) have carefully analysed their models for 
the sensitivity of the gravitational wave energy spectra to the 
neutron star radii and the adiabatic exponent of the polytropic 
equation of state. They investigated the possibility to deduce 
information from the heights and frequencies of the spectral 
peaks on these neutron star properties. Our energy spectra are 
in gross agreement with the principal features seen in their 
run 2 which simulates the coalescence of two neutron stars 
with 15 km radius. In particular, we recover the gradual drop 
of the spectrum below the point-mass inspiral value near the 
frequency at which the dynamical instability sets in, and the ex- 
istence of a primary peak associated with the main coalescence 
and of a secondary peak at higher frequency as a consequence 
of oscillations that occur during and after the merging. Also, 
the characteristic frequencies at which these peaks appear in 
our simulations match well with those of the Zhuge et al (1994) 
computation. The detailed structure of the energy spectrum, 
however, depends very sensitively on the total angular momen- 
tum of the system and therefore on the spins of the merging 
neutron stars. Moreover, the sharpness of the spectral structures 
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seen by Zhuge et al (1994) is barely reproduced by our results. 
This may be due to the inclusion of effects from gravitational 
radiation (back)reaction in our models, but can also be a con- 
sequence of the bulk viscosity of the nuclear medium which 
is caused by shifts of the /3-equilibrium and the correspond- 
ing production of neutrinos during compression and expansion 
phases. Note that since the disk is semi-transparent to neutri- 
nos and the central, massive object is not completely opaque, 
neutrinos can escape and are not in perfect thermodynamical 
equilibrium with the stellar gas. Neutrino production and emis- 
sion therefore give rise to dissipation. Both neutrino processes 
and gravitational radiation backreaction may be responsible for 
damping mass motions and smoothing fine structures. Finally, 
we point out that the shape of the emitted gravitational wave 
energy spectrum is strongly influenced by the duration of the 
intermediate phase of the stability of the merged object before 
its likely collapse into a black hole occurs. We find a drastic 
variation of the energy spectra with time and the spectra at 
about 3-5 ms after the start of our simulations (the expected 
time of the instability against gravitational collapse) are largely 
different from those obtained from integration until the end of 
our computations {t « 10 ms). 

Inspiral and final merger of binary neutron stars (each of 
the neutron stars could also be replaced by a black hole) is 
one of the most promising types of sources which might be 
detected by the gravitational wave interferometers that are cur- 
rently planned or under construction (see, e.g., Thorne, 1992). 
The maximum sensitivity of the LIGO experiment lies well 
below 1 KHz (namely around 100 Hz), which corresponds to 
an orbital distance of the neutron star binary of approximately 
35 km. Therefore most of the interesting signal is emitted dur- 
ing the inspiral phase when the orbit of the neutron stars is still 
larger than 35 km and can be described by using the point-mass 
approximation with a correction term to account for the finite 
extent of the neutron stars. When the binary separation is com- 
parable to the neutron star radius and the orbit gets close to the 
point of dynamical instability, hydrodynamic effects become 
important and the coalescence proceeds much different from 
the point-mass evolution. During the final stage of the merging 
the inspiraling binary neutron stars produce gravitational waves 
at dominant frequencies between 1 and 2 KHz. Although these 
frequencies lie outside the window of the largest sensitivity of 
the LIGO interferometers, the gravitational wave signals should 
nevertheless be detectable out to distances of several hundred 
Mpc. The emitted gravitational wave forms and spectra emitted 
during the merging of the binary are not only sensitive to the 
neutron star properties and the character of the equation of state 
of neutron star matter The hydrodynamics and mass motions 
depend on the angular momentum and, to a somewhat lesser 
extent, also on microphysical processes in the merging stars. To 
extract valuable information about the merging event from the 
measured signals therefore requires detailed numerical mod- 
els. In particular, steps away from the still basically Newtonian 
hydrodynamic models towards (fully) relativistic simulations 
need to be taken. 



In a forthcoming paper we shall present a detailed evaluation 
of our models for the neutrino emission and will study the 
implications for 7-ray burster models. Neutrino-antineutrino 
annihilation in the surroundings of the merging neutron stars 
during the last stages of the coalescence has been suggested 
as a possible mechanism to create relativistic e'^e fireballs 
which could be an efficient source of 7-rays (see, e.g., Cavallo 
& Rees, 1978; Paczyhski, 1990; Narayan et al, 1992; Meszaros 
& Rees, 1992; Mochkovitch et al, 1993). In future we intend 
to extend our simulations to cover coalescence of neutron stars 
with different masses, collisions of neutron stars that are initially 
not in orbit around each other, and also the merger of a neutron 
star with a black hole where we shall attempt to include the 
main effects of the radius of the last stable orbit. 

Movies in mpeg format of the dynamical evolution 
of all models are available in the WWW at http: 
/ /www . mpa-garching . mpg . de/~ mor/ nsgrb . html 

Acknowledgements. The calculations were performed at the Rechen- 
zentnim Garching on a Cray-YMP 4/64 and a Cray-EL98 4/256. It is a 
pleasure to thank Wolfgang Keil for transforming Lattimer & Swesty's 
FORTRAN equation of state into a usable table and M.R. would like to 
thank Sabine Schindler for her patience in our office. H.-Th. J. would 
like to thank W. Hillebrandt for inspiring discussions. H.-Th. J. was 
supported in part by the National Science Fovmdation under grant NSF 
AST 92-17969, by the National Aeronautics and Space Administra- 
tion under grant NASA NAG 5-2081, and by an Otto Hahn Postdoc- 
toral Scholarship of the Max-Planck-Society. M.R. acknowledges sup- 
port by the Deutsche Agentur fur Raumfahrtangelegenheiten (DARA), 
Forderungsvorhaben des Bundes, FKZ: 50 OR 92095. 

References 

Abramovici, A., et al, 1992, Science, 256, 325. 

Bildsten, L., Cutler, C, 1992, ApJ, 400, 175. 

Blanchet, L., Damour, T, Schafer, G., 1990, MNRAS, 242, 289. 

Bradaschia, C, et al, 1991, in Gravitation, Proc. Banff Summer Inst., 
Banff, Alberta; eds. R. Mann & P. Wesson; World Scientific, Sin- 
gapore. 

Bruenn, S.W., 1985, ApJS, 58, 771. 

Cavallo, G., Rees, M.J., 1978, MNRAS, 183, 359. 
Clark, J.PA., van den Heuvel, E.PJ., Sutantyo, W., 1979, A&A, 72, 
120. 

Colella, P, Woodward, PR., 1984, ICR 54, 174. 
Colella, P, Glaz, H.M., 1985, JCP 59, 264. 

Cooperstein, J., van den Horn, L.J., Baron, E., 1986, ApJ, 309, 653. 

Cooperstein, J., van den Horn, L.J., Baron, E., 1987, ApJ, 321, L129. 

Cutler, C, Flanagan, E.E., 1994, Phys. Rev., 49, 2658. 

Danzmann, K., etal, 1995, GEO600-A 600m Laser Interferometric 
Gravitational Wave Antenna, Proc. of the first Edoardo Amaldi 
Conference, Frascati, June 1994, to be published by World Scien- 
tific (Singapore). 

Danzmann, K.,etal, 1995, LISA -Laser Interferometer Space Antenna 
for gravitational wave experiments, Proc. of the seventh Marcel 
Grossmann Meeting, Stanford, July 1994, to be published by World 
Scientific (Singapore). 

Davies, M., Benz, W, Piran, T, Thielemann, F.K., 1994, ApJ, 431, 
742. 

Eastwood, J.W., Brownrigg, D.R.K., 1979, ICR 32, 24. 
Finn, L.S., Chemoff, D., 1993, Phys.Rev.D, 47, 2198. 



M. Ruffert et al.: Coalescing neutron stars - a step towards physical models 



31 



Finn, L.S., 1994, in Proc. of the Lanczos International Centenary 
Conference, eds. D. Brown et al, SIAM, Philadelphia. 

Fryxell, B.A., Muller, E., Arnett, D., 1989, Max-Planck-Institut fiir 
Astrophysik Report 449, Garching, Germany. 

Imshennik, VS., Popov, D.V., 1994, Astron. Lett., 20, 529. 

Iyer, B.R., Will, CM., 1993, Phys. Rev. Lett., 70, 113. 

Kochanek, C.S., 1992, ApJ, 398, 234. 

Kokkotas, K.D., Schafer, G., 1995, MNRAS, 275, 301. 

Lai, D., Rasio, F.A., Shapiro, S.L., 1994, ApJ, 420, 811. 

Lai, D., 1994, MNRAS, 270, 611. 

Lattimer, J.M., Swesty, RD., 1991, Nucl. Phys., A535, 331. 
Lincoln, C.W., Will, CM., 1990, Phys. Rev. D, 42, 1123. 
Meszaros, R, Rees, M.J., 1992, ApJ, 397, 570. 
Misner, C, Thome, K., Wheeler, J., Gravitation, 1973, Freeman, New 
York. 

Mochkovitch, R., Hemanz, M., Isem, J., Martin, X., 1993, Nature, 
361, 236. 

Nakamura, T, Oohara, K., 1991, Prog. Theor. Phys., 86, 73. 

Narayan, R., Paczynski, B., Piran, T, 1992, ApJ, 395, L83. 
Paczynski, B., 1990, in Proc. Los Alamos Workshop on Gamma-Ray 

Bursts, eds. H. Cheng et al, p. 67, Cambridge Univ. Press. 
Phinney, E.S., 1991, ApJ, 380, L17. 

Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1992, 

Numerical Recepies, Cambridge Univ. Press. 
Rasio, FA., Shapiro, S.L., 1992, ApJ, 401, 226. 
Rasio, F.A., Shapiro, S.L., 1994, ApJ, 432, 242. 
Reisenegger, A., Goldreich, P., ApJ, 426, 688. 
Ruffert, M., 1992, A&A, 265, 82. 

Shapiro, S.L., Teukolsky, S.A., Black Holes, White Dwarfs and 

Neutron Stars, Wiley, New York. 
Schinder, PJ., Schramm D.N., Wiita P.J., MargoUs S.H., Tubbs D.L., 

1987, ApJ, 313, 531. 
Shibata, M., Nakamura, T, Oohara, K., 1992, Prog. Theor. Phys., 88, 

1079. 

Shibata, M., Nakamura, T, Oohara, K., 1993, Prog. Theor. Phys., 89, 
809. 

Thome, K., 1987, in 300 Years of Gravitation, ed. S.Hawking & 

W.Israel, Cambridge Univ. Press. 
Thome, K., 1992, in Recent Advances in General Relativity, 

eds. A. Janis & J. Porter, Birkhauser, Boston. 
Tubbs, D.L., Schramm, D.N., 1975, ApJ, 201, 467. 
Zhuge, X., Centrella, J.M., McMillan, L.W., 1994, Phys.Rev, 50, 6247. 

A. Neutrino opacities 

In this appendix the formulae are collected to compute neutrino 
opacities /t = A^' (inverse mean free paths), optical depths, 
and diffusion timescales for all types of neutrinos and for the 
neutrino processes listed in Eq. (31)-Eq. (33). 

The transport cross section (cross section for momentum 
transfer) for neutrino-nucleon scatterings is given by 



(^sil^iN) = Cs,Ar (TO f — % ) 



(Al) 



where e is the neutrino energy, me is the electron rest 
mass, c the speed of hght, ao = 1.76 x 10~''^cm^, and 
Cs,n = (l + 5a^) /24 for neutrino-neutron scattering and 
Cs,p = [4(Cy — 1)^ + 5a^] /24 for neutrino-proton scattering 
with Cv = j + 2 sin^ ^w, sin^ ~ 0.23, and a « 1.25. 



We assume that the neutrino spectra can be represented by 
Fermi-Dirac distributions for the temperature T (taken equal to 
the gas temperature) and neutrino degeneracy r]^. = fi^jT 
being the neutrino chemical potential). 7]^^ is chosen to be 



riv^ = 



(A2) 



for heavy-lepton neutrinos and 

= €^ ■ [1 - exp(-T,.,o)] + nl^ ■ exp(-T,.,o) (A3) 
for electron neutrinos and 



= "C^ ■ [1 - exp(-Ti,^,o)] +ril^ ■ exp(-Ti,^,o) 



(A4) 



for electron antineutrinos. rf^ is the degeneracy parameter for 
Ve at chemical equiUbrium with the stellar medium, 



(A5) 



when rje is the degeneracy parameter of electrons (including 
electron rest mass contributions), rjp and rjn are the degeneracy 
parameters of p and n, respectively (without rest masses), and 
Q = 1 .2935 MeV is the rest-mass-energy difference between 
a neutron and a proton. To avoid divergent behavior at low 
densities we interpolate for i/g and i?e between the chemical 
equiUbrium values at high optical depths and fixed values r?°^ 
and ril in the limit of transparent matter. We use rj^ = rfl = 0. 
The interpolations between both limits are expressed in terms of 
the optical depths r^^ ,o and to^ ,o for Ve and i>e number transport 
(see below). Since the computation of optical depths already 
requires the knowledge of neutrino opacities, a first iteration 
step is performed where the deviation of r^^^ and 7]^^ from their 
chemical equiUbrium values is simply written as a function of 
decreasing density p. 

With defined neutrino spectra the spectraUy averaged scat- 
tering opacities are now computed as 

\meC^ ) T2+j{riyJ 

For j = one obtains the opacities for neutrino-number trans- 
port, for j = 1 the opacities for neutrino-energy transport. A is 
Avogadro's constant, T the gas temperature (measured in en- 
ergy units) and Tkirj) = dx x'' / [l+ exp(a; - r])] are the 
Fermi integrals for relativistic particles. Yi\t]\t means the num- 
ber fractions of free neutrons n or protons p, respectively, with 
the double index "NN" indicating that Pauli blocking effects 
for degenerate nucleons are taken into account by following 
Bruenn (1985), who suggested the replacements 



_ JFat for nondegenerate TV ; 

- <^YN/jm for degenerate TV . ^^'^ 



Simple interpolation of both limits yields the combined formula 



f JVJV 



1 + |max(?;jv,0) 



(A8) 
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For completely dissociated matter the nucleon fractions are 
Y„ = l-Y,mdYp = Y,. 

The cross sections for absorption onto n in hot neutron 
star matter is given by 



(TAVen) = (To . 

and for absorption onto p it is 



l+3a^ 

(TAVeP) = -. Cro 



[1 -/e-(e)] 



[1 -/e-(e)] 



(A9) 



(Tubbs & Schramm, 1975; Bruenn, 1985), where the factors in 
brackets are supposed to account for Pauli blocking effects in 
the electron and positron phase spaces, respectively. Performing 
the spectral averaging one finds the opacities for absorption 
processes as 



, , l + 3a2 

Kajil^eTl) = aoApYnp 



(l-/e-(e)) 



(All) 



and 



Ka,j(j^eP) = — — aoApYpn 

•^2+j(f?sJ 



(1 - fe^e)) 



(A12) 



Again, the doubly indexed number fractions Ynp and indi- 
cate the inclusion of Fermion blocking effects in the nucleon 
phase spaces (Bruenn 1985) and are defined as 



2Y, - 1 



exp(??p - ??n) - 1 

and 

Ypn exp(7^j, 7^^) ■ Yjip 



(A13) 



(A14) 



when Yn = I — Ye and = are used. For the phase space 
blocking of the leptons we employ the approximate expressions 



(l-/e-(e)) = {l+exp 
and 

(l-fAe)) = |l+exp 



f ^5iv,J , All 



-1 



(A15) 



(A16) 



From the spectral averages of scattering and absorption 
opacities total transport opacities are computed for Ug, 

KljiVe) = KsjiVeTl) + Ksji^eP) + KnjiVen) , (A17) 
for l>e, 

Kt,j(z>e) = /^s,j(^'e^^) + Ks,j(z>eP) + Ka,j(l^eP) , (A18) 



and for the heavy-lepton neutrinos u^. 



(A19) 



These allow us to compute optical depths along paths [si, S2] 
according to 



T,,i,ji[s 1,82]) 



/S2 
-1 



ds Kfji^i) 



(A20) 



(A 10) which then yields the diffusion timescale along the path as 



tf^J ([SUSI]) « ^^^'^ T,„j([suS2]) . (A21) 

The numerical factor in Eq. (A21) is of order unity but its ex- 
act value cannot be fixed from theoretical arguments. It must 
therefore be adjusted in course of a calibration of the neu- 
trino leakage scheme with more precise transport methods. In 
three-dimensional situations the optical depths and diffusion 
timescales for all cells of the grid into all directions (±x, ±y, 
±z) are evaluated. The diffusion timescale attributed to a par- 
ticular cell is then chosen as the minimum of these diffusion 
times. 



B. Neutrino emission rates 

Given the neutrino diffusion timescales the effective neutrino 
emissivities of all cells of our computational grid can be deter- 
mined for the processes of Eq. (27)-Eq. (30). 

The emission rate of Ug (cm^^s ') by the /3-process of 
Eq. (27) in hot neutron star matter is given by 

Rpii^e) = '^0'^ ApYp„ ■ Eg- ■ (1 - (Bl) 

o \TYlgC ) 

and the corresponding emission rate of in process Eq. (28) 
is 



1 + 3a^ (tqc 
8 (mec2)2 



ApYn 



(B2) 



(Tubbs & Schramm, 1975; Bruenn, 1985). Y^p and Yp„ are de- 
finedinEq. (A13)andEq. (A14), respectively, and the blocking 
in the neutrino phase spaces is approximately taken into account 
by 



(l-/.e(e))/3 
and 

{l-IoMp = |l + exp 



(B3) 



where the Fermi functions are evaluated with the average ener- 
gies of captured electrons or positrons (approximately equal to 
the energies of produced or z>e). In Eq. (Bl)andEq. (B2)and 
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in the rates given further below energy moments of the electron 
and positron distributions occur, being defined as 



and 



(/ic)3 

Stt 
(/ic)3 



(hey 



(B5) 
(B6) 

(B7) 



The emission of or i/g by electron-positron pair annihi- 
lation (Eq. (29)) is given by 



36 {med-y 
(1 -/..(e))e. (1 -/^.(e)).e 



(B8) 



(Cooperstein et al, 1986, 1987) with the Ve and phase space 
blocking factors being evaluated with the average energies of 
neutrinos produced by e'^e -annihilation: 

(l-/.,(e))ee = 

[ (I ^4(r?e) . 1 ^4(-r?e) \\\~' 

The weak interaction constants are (C\ + €2)^ ,5 = {Cy — 
Ca)^ + {Cv + Ca)^ with Ca = j- For the production of j/^, Vr, 
D^, and Pt- the corresponding rate is 



(Ci + C2)v^o^ apc 
9 (mec2)2 

((l-/..(e))ee)' 



(BIO) 



with (Ci + C2),,^, = (Cv - C^)2 + {Cv + Ca- 2)2. 

The rate of creation of Ve or Dg by the decay of transversal 
plasmons (Eq. (30)) (longitudinal plasmons yield a negligible 
contribution, see Schinder et al. 1987) can be written with suf- 
ficient accuracy as 



ia* 

7'^e-^(l +7) ■ (1 - US^))^ (1 - hS^% (Bll) 
and the corresponding rate for producing Vfj_, v^,, Vr, and Vr is 

7'e-^(l+7)- ((l-/..(e))-y)' • (B12) 

a* is the fine-structure constant, a* = 1/137.036, and 7 « 

70 (7r2 + 3r?2) with 70 being related to the plasma frequency 

by 70 = hfla/mefP' = 2 \/a*~l(iTr) = 5.565 • 10^^. The average 
blocking factor for Pi production can be estimated by evaluating 



the Fermi function f,,. (e) with the average energy of neutrinos 
created by the plasmon decays: 

(1 - UM-r = {1 + exp - (^1 + - r,,)j | . (B13) 

The neutrino-energy emission rates by the /3-processes, 
e'^e~-pair annihilation, and plasmon decay are 



Eg* 



QeeW) = ReeiVi) ' ^ 



and 



1 



Q^(Pi) = RJvi)--T 2 + 



1+7 



(B14) 
(B15) 
(B16) 

(B17) 



for neutrino species Vi. 

Making use of the neutrino number densities. 



At: 



9ui 



{hcf 

and energy densities. 



9v, 77-3 T'-Tiinv,)^ 
{hey 



(B18) 



(B19) 



where the numerical multiplicity factor is g^^ = g^^ = I and 
g^^ = 4, we can define neutrino-number emission timescales 
as 

te)"' = — [Rp{Vi) + Ree{Vi) + R^{Vi)] = (B20) 

and neutrino-energy emission timescales t^°^\ as 

{t'^^^.y' = — [Qp{Vi)+Qee{l^i) + Qj{l^i)\ = ^ (B21) 

(of course, the /^-processes only apply for Vg and Vg). Effective 
emission rates of neutrino number and energy can now be de- 
fined as functions of the relative size of the shortest diffusion 
timescale of a grid cell and the neutrino emission timescale 
t^°f - of this cell: 



R''\vi) 



R{vi) 



^ + ,0 y-v, ,0 ) 
Qivi) 



-1 ' 



l + t 



diff ( ^loss ( 



(B22) 



(B23) 



In the (nearly) transparent, low-density regime the optical 
depth to neutrinos vanishes and the diffusion timescale is short 



34 



M. Ruffert et al.: Coalescing neutron stars - a step towards physical models 



compared to the timescale of direct neutrino loss. Therefore 
R'^^i'i) R{vi) and Q^'^^Vi) Q{vi). On the other hand, in 
the opaque region the diffusion timescale becomes very long 
(and due to the increase of the rates with density and tempera- 
ture the direct loss timescale becomes extremely short, too). In 
this case R^^ivd n^, /tf^g and Q^'^Vi) e^, /tf^i, i.e. the 
loss of neutrinos is then determined by the much slower diffu- 
sion process that depletes the equilibrium neutrino distributions 
only on a long time. 

Finally, we define the total electron-lepton number loss rate 
of the stellar gas as 

Ry^ = - R^^ive) + R^^\i>e) (B24) 
and the total energy loss rate as 

Q; = - [Q'^^iue) + Q'\ve) + Q'\v^)\ . (B25) 

Eq. (B24) and Eq. (B25) are the source terms to be used in 
the continuity equation for electron-lepton number (Eq. 14: 
S-L = Ry ) and in the gas-energy equation (Eq. 3: Se = Qj), 
respectively. 

The average energy of neutrinos of species emitted from 
a grid cell is given by 



= m^) = Ri^) — 



and the average energy of neutrinos Vi emitted from the star is 
determined as the ratio of the rates Q^^ivi) and R^{vi) sunmied 
individually over the whole grid. 
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